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Abstract 

The time fractional diffusion equation is obtained from the standard 
diffusion equation by replacing the first-order time derivative with a 
fractional derivative of order /3 £ (0, 1) . The fundamental solution for 
the Cauchy problem is interpreted as a probability density of a self- 
similar non-Markovian stochastic process related to a phenomenon of sub- 
diffusion (the variance grows in time sub-linearly). A further generaliza- 
tion is obtained by considering a continuous or discrete distribution of 
fractional time derivatives of order less than one. Then the fundamental 
solution is still a probability density of a non-Markovian process that, 
however, is no longer self-similar but exhibits a corresponding distribution 
of time-scales. 
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1 Introduction 



The main physical purpose for adopting and investigating diffusion equations 
of fractional order is to describe phenomena of anomalous diffusion usually met 
in transport processes through complex and/or disordered systems including 
fractal media. In this respect, in recent years interesting reviews, see e.g. [2^1 
1311 140j. have appeared, to which (and references therein) we refer the interested 
reader. All the related models of random walk turn out to be beyond the classical 
Brownian motion, which is known to provide the microscopic foundation of the 
standard diffusion, see e.g. P^[57] . 

The diffusion-like equations containing fractional derivatives in time and/or 
in space are usually adopted to model phenomena of anomalous transport in 
physics, so a detailed study of their solutions is required. 

Our attention in this paper will be focused on the time fractional diffusion 
equations of a single or distributed order less than 1, which are known to be 
models for sub-diffusive processes. Specifically, we have worked out how to 
express their fundamental solutions in terms of series expansions obtained from 
their representations with Mellin-Barnes integrals. 

In Section 2 we shall recall the main results for the fundamental solution of 
the time fractional diffusion equation of a single order, which are obtained by 
applying two different strategies in inverting its Fourier-Laplace transform. Both 
techniques yield the same power series representations of the required solution: 
it turns out to be self similar (through a definite space-time scaling relationship) , 
and expressed in terms of a special function of the Wright type. Then, in Section 
3, we shall apply the second strategy for obtaining the fundamental solution in 
the general case of a distributed order. We have provided a representation in 
terms of a Laplace-type integral of a Wright function, that can be expanded 
in a series containing powers of the space variable and certain functions of 
time, responsible for the time-scale distribution. Finally, in Section 4, the 
main conclusions are drawn. For convenience and self-consistency we provide 
an Appendix devoted to our notations of fractional calculus. 

2 The time fractional diffusion equation of 
single order 

The standard diffusion equation in re-scaled non-dimensional variables has the 
form 

d (9^ 

— u{x,t):^—u{x,f), xeR, teRc[, (2.1) 

with u{x,t) as the field variable. We assume the initial condition 

«(a;,0+) =uo(x), (2.2) 

where uq{x) denotes a given ordinary or generalized function defined on R, 
that is Fourier transformable in ordinary or generalized sense, respectively. We 
assume to work in a suitable space of generalized functions where it is possible 
to deal freely with delta functions, integral transforms of Fourier, Laplace and 
Mellin type, and fractional integrals and derivatives. 
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It is well known that the fundamental solution (or Green Junction) of Eq. (2.1) 
i.e. the solution subjected to the initial condition u{x,Q) = uq{x) = 5{x), and 
to the decay to zero conditions for — > oo, is the Gaussian probability density 
function [pdf) 

u(x,t) = -^t-V2e--V(40, (2.3) 
that evolves in time with second moment growing linearly with time, 

/+ca 
x'^u{x,t)dx = 2t, (2.4) 
-oo 

consistently with a law of normal dijfusioi^. We note the scaling property of 
the Green function, expressed by the equation 

u{x, t) = t"i/2 u{x/t^''^) , with U{x) := u{x, 1) . (2.5) 

The function U{x) depending on the single variable x turns out to be an even 
function U{x) = U{\x\) and is called the reduced Green Junction. The variable 
X := x/t^^^ is known as the similarity variable. It is known that the Cauchy 
problem {(2.1) — (2.2)} is equivalent to the intcgro-diffcrential equation 



w(x, t) — uo{x) 
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dr , (2.6) 



where the initial condition is incorporated. Now, by using the tools of the 
fractional calculus we can generalize the above Cauchy problem in order to 
obtain the so-called time fractional diffusion equation in the two distinct (but 
mathematically equivalent) forms available in the literature, where the initial 
condition is understood as (2.2). For the essentials of fractional calculus we 
refer the interested reader to the Appendix. 

If /? denotes a real number such that < /? < 1 the two forms are as follows: 
d 5^ 

—u{x,t)^tD^-f^ — u{x,t), a;GR,teR([; u{x,Q+) ^ uq{x) , (2.7) 

where tD^~^ = tD^ tJ^ denotes the Riemann-Liouville (R-L) time-fractional 
derivative of order 1 — /?, see (A. 4) with to = 1 , and 

tD^uix,t) = —u{x,t), a;eR,<eR+; u{x,0+) = uoix) , (2.8) 
ox'^ 

where tD^ = tJ^~^tD^ denotes the time-fractional derivative of order f3 
intended in the Gaputo sense, see (A. 5) with to = 1. In analogy with the 
standard diffusion equation we can provide an integro-differential form that 
incorporates the initial condition (2.2): for this purpose we replace in (2.6) the 



^The centred second moment provides the variance usually denoted by u^{t). It is a 
measure for the spatial spread of u{x,t) with time of a random walking particle starting at 
the origin a; = 0, pertinent to the solution of the diffusion equation (2.1) with initial condition 
u{x, 0) = S{x). The asymptotic behaviour of the variance as t — > oo is relevant to distinguish 
normal diffusion {a-^{t)/t — > c > 0) from anomalous processes of sub-diffusion (a^{t)/t — > 0) 
and of super-diffusion icr'^{t)/t — > -l-oo). 
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ordinary integral with the Riemann-Liouville fractional integral of order (3, tJ^, 
namely 

r 5^ 

u{x,t) = uo{x) + tJ^ ^n{x,t . (2.9) 



Then, the above equations read explicitly: 



— u(a;,r) 



dr 



u{x, r) 



dr 



it-r) 



u{x, 0^) = uo{x) . 



(2.7') 



dr 



u{x, t) = 1*0(2;) 



— u{x, t) , u{x, 0+) = uoix) , (2.8') 



1 



m Jo 



t r a2 



dx"^ 



u(x, t) 



dr 



(t-r) 



(2.9') 



The two Cauchy problems (2.7), (2.8) and the integro-differential equation (2.9) 
are equivalenci: for example, we derive (2.7) from (2.9) simply differentiating 
both sides of (2.9), whereas we derive (2.9) from (2.8) by fractional integration 
of order p. In fact, in view of the semigroup property (A. 2) of the fractional 
integral, we note that 

tJ^ tD^uix, t) = tJ^ tJ^-l" tD\{x, t)=tJ^ tD\ix, t)^u{x, t) - uo{x). (2.10) 

In the limit /3 = 1 we recover the well-known diffusion equation (2.1). 

Let us consider from now on the Eq. (2.8) with uo{x) = 5{x): the fundamental 
solution can be obtained by applying in sequence the Fourier and Laplace 
transforms to it. We write, for generic functions v(x) and u'(i), these transforms 
as follows: 



T {v{x)] k\ = v{Hi) := J 
C{w{t);s} = w{s) 



+ 00 



-st, 



v{x) dx , K G R , 
!{t) dt, see. 



(2.11) 



Then, in the Fourier-Laplace domain our Cauchy problem [(2.8) with u(x, 0^) = 
5{x)], after applying formula (A. 6) for the Laplace transform of the fractional 
derivative and observing (5(k) = 1, see e.g. |10j . appears in the form 'u{k, s) — 



-K^ u{k, s) , from which we obtain 



u(k, s) 



0</3<l, 5R(s)>0, kgR. 



(2.12) 



To determine the Green function u(x, t) in the space-time domain we can follow 
two alternative strategies related to the order in carrying out the inversions in 
(2.12). 



^The integro-differential equation (2.9) was investigated via Mellin transforms by Schneider 
& Wyss |35 | in tlieir pioneering 1989 paper. The time fractional diffusion equation in the form 
(2.8) with the Caputo derivative has been preferred and investigated by several authors. From 
the earlier contributors let us quote Caputo himself [3], Mainardi, see e.g. 1201 1211 1221 and 
Gorenflo et al. |16l 117) . In particular, Mainardi has expressed the fundamental solution in 
terms of a special function (of Wright type) of which he has studied the analytical properties 
and provided plots also for 1 < /3 < 2, see also | 12lll3ll25l and references therein. For the form 
(2.7) with the R-L derivative earlier contributors include the group of Prof. Nonnenmacher, 
see e.g. |28| . and Saichev & Zaslavsky I33| . 
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(51) : invert the Fourier transform getting u(x, s) and then invert the remaining 
Laplace transform; 

(52) : invert the Laplace transform getting u{k, t) and then invert the remaining 
Fourier transform. 

Strategy (SI): The strategy (SI) has been apphed by Mainardi [201 HD [12] to 
obtain the Green function in the form 

u{x,t) ^t-^''^U {\x\/tf^/^^ , -cx)<x<+cx), i>0, (2.13) 

where the variable X xjt^l"^ acts as similarity variable and the function 
U{x) :— u{x,l) denotes the reduced Green function. Restricting from now on 
our attention to a; > 0, the solution turns out as 

1 -j^ / x)^ 

C/(x) = -M|(x)= 2T.klT[-pk/7+ (1-^/2)] = 

^ E r[(/3(fc + l)/2] sin[(7r/3(fc + l)/2] , 

fc=0 

(2.14) 

where M£{x) is an an entire transcendental function (of order 1/(1 — of 
the Wright type, see also P^ITB] and [5^. 

Strategy (S2): The strategy (S2) has been followed by Gorenflo et al. pjj and by 
Mainardi et al. 24J to obtain the Green functions of the more general space-time 
fractional diffusion equations, and requires to invert the Fourier transform by 
using the machinery of the Mellin convolution and the Mellin-Barnes integrals. 
Restricting ourselves here to recall the final results, the reduced Green function 
for the time fractional diffusion equation now appears, for a; > 0, in the form: 

U{x) = - cos {nx) E, d.= -— J^_^^ ril-Ps/2) ^ ' ' 

(2.15) 

with < 7 < 1, where Ep denotes the Mittag-Lefflcr function, see e.g. [51 123j. 
By evaluating the Mellin-Barnes integrals using the residue theorem, we arrive 
at the same power series (2.19). 

Both strategies allow us to prove that the Green function is non-negative and 
normalized, so it can be interpreted as a spatial probability density evolving in 
time with the similarity law (2.13). Although the two strategies are equivalent 
for yielding the required result, the second one appears more general and so more 
suitable for treating the more complex case of fractional diffusion of distributed 
order, see the next Section. 

Of particular interest is the evolution of the second moment as it can be derived 
from Eq. (2.12) noting that 

_ ^2 _ 2 t/^ 

When < /? < 1 the sub-linear growth in time is consistent with an anomalous 
process of sub-diffusion. 
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3 The time- fractional diffusion equation of dis- 
tributed order 



The fractional diffusion equation (2.8) can be generalized by using the notion of 
fractional derivative of distributed order in tim^. We now consider the so-called 
time-fractional diffusion equation of distributed order 

^ b{/3) [ tD^ u{x, t)] d/3 u{x, t) , b{(3) > 0, ^ b{(3) d/3 = 1 , (3.1) 

with X € R, i > 0, subjected to the initial condition u(a:, 0+) — 5{x). Clearly, 
some special conditions of regularity and boundary behaviour will be required 
for the weight function fe(/3), that we call the order-density. 

Equations of type (3.1) have recently been discussed in [6l El [3 [36] and in [30] . 
As usual, we have considered the initial condition ^(a;, 0"*") = 5{x) in order to 
keep the probability meaning. Indeed, already in the paper [B], it was shown that 
the Green function is non-negative and normalized, so allowing interpretation 
as a density of the probability at time t of a diffusing particle to be in the 
point X. The main interest of the authors in [T] [3 [M] was devoted to the 
second moment of the Green function (the displacement variance) in order to 
show the subdiffusive character of the related stochastic process by analyzing 
some interesting cases of the order-density function &(/3). For a thorough general 
study of fractional pseudo-differential equations of distributed order let us cite 
the paper by Umarov and Gorenflo [39j . For a relationship with the Continuous 
Random Walk models we may refer to the paper by Gorenflo and Mainardi [15] . 



In this paper, extending the approach by Naber [3D], we provide a general 
representation of the fundamental solution corresponding to a generic order- 
density &(/3). By applying in sequence the Fourier and Laplace transforms to 
Eq. (3.1) in analogy with the single-order case, see Eq. (2.12), we obtain 



from which 



where 



b{(i) d(i 
u{k, s) ~ 



u{k,s)- I b{(3)s^-^ dl3 ^ -k^ u{k,s) , 
B{s)/s 



B{s) -f k2 



B{s) 



(s) > , K e R, 



&(/?) sl^dp. 



(3.2) 
(3.3) 



Before trying to get the solution in the space-time domain it is worth to outline 
the expression of its second moment as it can be derived from Eq. (3.2) using 
(2.16). We have, for k, near zero, 

= K' " ^ ^ ■ ■ ■) ' ^'^"^ ^ ^ °' ^ ^^-^^ 

^We find an earlier idea of fractional derivative of distributed order in time in the 1969 
book by Caputo [3], that was later developed by Caputo himself, see [4|[5] and by Bagley & 
Torvik, see [l]. 
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Then, from (3.4) we are allowed to derive the asymptotic behaviours of /i2(t) for 
t ~-^ 0^ and t — *■ +00 from the asymptotic behaviours of B(s) for s — > 00 and 
s ^ 0, respectively, in virtue of the Tauberian theorems. The expected sub- 
linear growth with time is shown in the following special cases of &(/3) treated 
in O [7] . The first case is slow diffusion (power-law growth) where 

h{(3) = hx5{(3 - (3i) + h25{l3 ~ P2), < /3i < ^2 < 1, 61 > 0, 62 > 0, 61 + &2 = 1- 

In fact 

tp\ t 



~i ^ 2 I 62r(/32 + 1) 

6ir(/3i + 1) 



In [3], see Eq. (16), the authors were able to provide the analytical expression 

of /i2(i) in terms of a 2-parameter Mittag-Lcfflcr function. 

The second case is super-slow diffusion (logarithmic growth) where 

6(/3) = 1, 0</3<l. 

In fact 



logs 



2nog(lA), 



^(s)=2— so ^i2{t)-{ (3.6) 



s(s-l) 



2 log(t), 



In [5], see Eqs. (23)- (26), the authors were able to provide the analytical 
expression of /i2(0 in terms of an exponential integral function. 
Inverting the Laplace transform, in virtue of a theorem by Titchmarsh we obtain 
the remaining Fourier transform as 

u{K,t)^-^J^ e-''*Im|S(re*'^)} dr, (3.7) 

that requires — Im {B{s)/[s{B{s) + k^]} along the ray s = re"^ with r > (the 
branch cut of the functions s^ and s^^^). By writing 

( p = pir)^\B{re^^)\, 

B re*^ =pcos(7r7)-|-zpsin(^7), <^ 1 (3.8) 

V J I 7 = 7(r) -arg (re'^^jj, 

we get after simple calculations 

u{K,t) ^ K{K,r)dr, (3.9) 



where 

^ TT + 2k2 p cos(7r7) + p2 ' 

Since u{x, t) is symmetric in x, Fourier inversion yields 

1 /■+°° f e~^* 1 
u{x,t) ^ — / cos(kx) < / K{K,r)dr) dn . (3.11) 

TT In \ . n r 
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To calculate this Fourier integral we use the Mellin transform. Let 

M{fiO;s}^ ris)^ fiOC^dt 7i<5R(s)<72, (3.12) 



be the Mellin transform of a sufficiently well-behaved function /(^) , and let 

M-'{nsy,^} = f{0^,- ns)r'ds, (3.13) 

be the inverse Mellin transform, where > , 7 = 5R(s) , 71 < 7 < 72 ■ 

Denoting by <^ the juxtaposition of a function /(^) with its Mellin transform 
/*(s), the Mellin convolution theorem implies 



00 -j^ 



Jo V 

Then, following [25, pp. 160-161, we identify the Fourier integral in (3.11) as a 
Mellin convolution in k , that is u{x, t) — f{K, t) (g) g{K, x), if we set (see (3.14) 
with ^ = l/x, rj = k) 

poo -rt 

f{K,t):^ KiK,r)dr n ris,t), (3.15) 

Jo ^ 

g{K,x) -.^ ^—cosf-] <^ ^ ^sin(^) ■.^g*{s,x), (3.16) 

TT X K V'^/ ^ \ 2 / 

with < 5R(s) < 1. The next step thus consists in computing the Mellin 
transform f*{s,t) of the function f{K,t) and then inverting the product 
f*{s,t) g*{s,x) using (3.16) in the Mellin inversion formula, namely 

1 1 ra+icc 

u{x,t)= — — / /*(s,t)r(l - s)sin(7rs/2)a;"ds 

Trx2mj^_j^^ ('01 7^ 

,. r(i-.) '> 

- x2mj^_,^ ^ ^'''^r(./2)r(l-,s/2)'^ 
The required Mellin transform f*{s;t) is 



By the variable change — > p/i the term in braces becomes 



,./2+l 1 ^00 .,(./2+l)-l 



2p TT p? + 2/^cos(7r7) + 1 

(3.19) 

p'^/^ r r(s/2 + 1) r[l - (s/2 + 1) 



2 i r(7s/2)r(l-7s/2) 

where we use a formula from the Handbook by Marichev, see [27] p. 156, Eq. 
15 (1), under the condition < 3fi(s/2 + 1) < 2, I7I < 1. As a consequence of 
(3.18)-(3.19) we get 

r( A r ""'^ P'''^ f r(./2 + i)r[i-(./2 + i)] ] 

^^'^'^ = ~L r(7./2)r(i-7./2) ^'-''^ 
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Now, using Eqs (3.17) and (3.20) we can finally write the solution as 

1 r°° p~rt 

u{x, t)^— F[p^'^x) dr , (3.21) 

where F{p^^'^x) is expressed in terms of Mcllin-Barnes integrals: 

27rzy_ r(7./2)r(l-7./2)^^ 
— / T{1- s)sm{TT^s/2){p^'^xy ds , 



(3.22) 



with p = p{r), 7 = 7(r). We remind that p{r) and 7(r) are related to the 
order-density b{(3) according to Eqs. (3.3) and (3.8). By evaluating the Mellin- 
Barnes integrals via the residue theorem we arrive at the series representations 
in powers of (p^^^x), 

1/2 ^ ^ 1/2 ^ ^ 

' ^ ^£^^A:!r(7fc/2 + 7/2)r(-7A:/2 + 1-7/2) 

°° ('„ 1/2 \k (0.23) 

P""^T. \i sin(^7(fc + l)/2). 

fc=0 

This function can be interpreted again in terms of generalized Wright or H- 
Fox functions, as outlined in [26]. From Eqs. (3.21) and (3.23), interchanging 
integration and summation, we get the series representation of the fundamental 
solution: 

i^x^ 



k=0 

where, with p ~ p{r), 7 — 7(7"), 

r°° p-rt 

ipkit)= sin[^7(A: + l)/2]p(^-+i)/2(ir. (3.25) 



In order to check the consistency of this general analysis we find it instructive to 
derive from Eqs. (3.24)-(3.25) the results of Section 2 concerning the fractional 
subdiffusion of a single order. We denote this (fixed) order by ly in distinction 
from the variable /3 in the distributed order case. This means to consider in Eq. 
(3.1) the particular case 

b{(3) = S{P - ly) , 0<J/<1, (3.26) 

so that B{s) = s" and Eq. (3.8) yields 

p = p{r) = r"^ , 7 = const = v . (3.27) 

In the series representation of the fundamental solution (3.24)-(3.25) the 
functions (pk (t) now turn out as 



ifikit)^ siii[TTiy{k + l)/2] r^C" + '^y^ dr 

sin[7ri^(fc + l)/2] 



rli^jk + l)/2] 

^i>(fc+l)/2 



(3.28) 
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As a consequence, the solution reads 



u{x, t) = -t-^l'' ■ - J2 ^ ' r[Kfc + l)/2] sin[7rz.(fc + l)/2] = 

2 TT -^-^ fc! 

2 2 Vi^/sy 

(3.29) 

in agreement with Eqs. (2.13)-(2.14). Of course, only in this special case it is 
possible to single out a common time factor {t^'^^^) from all the functions (pk{t) 
and get a self-similar solution. In general the set of functions ipk (t) gives rise to 
a distribution of different time scales related to the order density b{j3). 

4 Conclusions 

After outlining the basic theory of the Cauchy problem for the spatially 
one-dimensional and symmetric time fractional diffusion equation (with its 
main equivalent formulations), we have paid special attention to transform 
methods for finding its fundamental solution or (exploiting self-similarity) the 
corresponding reduced Green function. We have stressed the importance of the 
transforms of Fourier, Laplace and Mellin and of the functions of Mittag-Leffler 
and Wright type, avoiding however the cumbersome _ff-Fox function notations. 

A natural first step for construction of the fundamental solution consists in 
applying in either succession the transforms of Fourier in space and Laplace 
in time to the Cauchy problem. This yields in the Fourier-Laplace domain 
the solution in explicit form, but for the space-time domain we must invert 
both transforms in sequence for which there are two choices, both leading to 
the same power series in the spatial variable with time-dependent coefficients. 
The strategy, called by us the "second", of first doing Laplace inversion and 
then the Fourier inversion yields the reduced Green function as a Mellin-Barnes 
integral form which, by the calculus of residues, the power series is obtained. 
This strategy can be adapted to the treatment of the more general case of 
the distributed order time fractional diffusion equation. Now the fundamental 
solution can be expressed as an integral over a Mellin-Barnes integral containing 
two parameters having the form of functionals of the order-density. Finally, 
again for the fundamental solution a power series comes out whose coefficients, 
however, are time-dependent functionals of the order-density. But, if there is 
more than one time derivative-order present, self-similarity is lost. 

Appendix. The two fractional derivatives 

For a sufficiently well-behaved function f{t) {t E R+) we may define the 
fractional derivative of order P {m — 1 < (3 < m , me N), see e.g. [T4|, [32] 
in two different senses, that we refer here as to Riemann-Liouville derivative 
and Caputo derivative, respectively. Both derivatives are related to the so- 
called Riemann-Liouville fractional integral of order a > 0, see defined 
as 

t^"/(i):-f^ f\t-Tr-'f{T)dT, a>0. (Al) 
r(a) Jo 
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We note the convention tJ^ = I (Identity) and the semigroup property 

t.r ^ ,r ^ tr+p , a>o, /5>o. (A2) 

The fractional derivative of order /3 > in the Riemann-Liouville sense is defined 
as the operator which is the left inverse of the Riemann-Liouville integral 
of order (3 (in analogy with the ordinary derivative) , that is 

tD^tJ^=I, /3>0. (A3) 



If m denotes the positive integer such that m ~ 1 < f3 < m , we recognize from 
Eqs. (A.2) and (A.3) tD^ f{t) := tD"" tJ"^-^ f{t) , hence 




fir) dr 



V{m- (3) Jo (t-T)/3+i-'" 
fit), f3 



m — 1 < /3 < TO, 

(^.4) 



For completion we define tZ?^ — I . On the other hand, the fractional derivative 
of order /3 > in the Caputo sense is defined as the operator such that 
tD^fit) tJ'»-^i?"V(^), hence 

' ■' ^ ' - l< (3 <m, 



,^e/(t)= r(m-/5)7o ' (A5) 

We point out the major utility of the Caputo fractional derivative in treating 
initial-value problems for physical and engineering applications where initial 
conditions are usually expressed in terms of integer-order derivatives. This can 
be easily seen using the Laplace transformation, according to which 

m — 1 

C{tD^Jit)-s)=sPJis)-Y,sP-^-'^ f(^\Q+), m-l</3<TO, (A6) 

poo 

where /(s) = £{/(t);s} = / e-^*f{t)dt, s e C and /(''■)(0+) := 

Jo 

lim f^'^Ut). The corresponding rule for the Riemann-Liouville derivative is 
more cumbersome: for m — 1 < P < m it reads 

m — 1 

£ {tD^ fit); s} = lis) - ^ [tD'^ /(0+) /"-^-''^ , (A.7) 

fc=0 

where, in analogy with (A. 6), the limit for t ^ 0+ is understood to be taken 
after the operations of fractional integration and derivation. As soon as all 
the limiting values /'•'"•' (0+) are finite and to — 1 < /3 < m, the formula (A. 7) 
simplifies into 

C{tDP fit);s]=sP fis). (A8) 

In the special case /(fc)(0+) for fc = 0, 1,TO- 1, we recover the identity 
between the two fractional derivatives. The Laplace transform rule (A. 6) 
was practically the starting point of Caputo |2t in defining his generalized 
derivative in the late sixties. For further reading on the theory and applications 
of fractional calculus we recommend the recent treatise by Kilbas, Srivastava & 
TrujiUo Hg. 
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